A microscopic mean-field theory of the jamming transition 
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Dense particle packings acquire rigidity through a nonequilibrium jamming transition commonly 
observed in materials from emulsions to sandpiles. We describe athermal packings and their observed 
geometric phase transitions using equilibrium statistical mechanics and develop a fully microscopic, 
mean-field theory of the jamming transition for soft repulsive spherical particles. We derive analyti- 
cally some of the scaling laws and exponents characterizing the transition and obtain new predictions 
for microscopic correlation functions of jammed states that are amenable to experimental verifica- 
tions, and whose accuracy we confirm using computer simulations. 



About 50 years ago, Bernal [H used dense disordered 
sphere packings as model systems to understand the liq- 
uid state, at a time where the statistical mechanics of liq- 
uids was still in its infancy. Today, the idea that jammed 
materials share deep similarities with dense liquids and 
glasses remains popular Q. However, while liquid state 
theory grew as a cornerstone of theoretical physics [1| , no 
equivalent theory is available for jammed matter, because 
this is a nonequilibrium, amorphous, athermal state of 
matter — a theoretical challenge overlooked by Bernal. 
Thus, despite intense research activity with a large 
body of numerical and experimental observations 0,13: it 
is not yet clear what is the appropriate theoretical frame- 
work to understand dense athermal packings and the in- 
triguing phase transitions they undergo, although foams, 
pastes and emulsions are familiar materials. 

We address the purely geometrical packing problem of 
soft spheres, and suggest to study first their statistical 
mechanics at finite temperatures, T, before taking the 
T — > limit where jamming occurs. A similar approach 
is frequently used in combinatorial optimization prob- 
lems [9|] , because powerful statistical mechanics tools can 
then be used in a context where they are not a priori 
relevant [T^]. We investigate the statistical mechanics of 
the system at T > using mean- field theory [llj , and 
develop a fully microscopic theoretical scheme to predict 
the structure of non-equilibrium configurations of soft re- 
pulsive spheres at zero [13] and finite [l^] temperatures. 

Our microscopic approach is thus markedly different 
from recent theoretical works [1, Q , which are based on 
phenomenological and scaling considerations. Similarly 
to Landau- Ginzburg theory of phase transitions, our aim 
is to derive, from first principles, the correct qualitative 
description of the transition, and accurate quantitative 
predictions for several observables. However, as a mean- 
field theory, our approach does not describe well all fluc- 
tuations near the transition, and the associated scaling 
laws 0, 0]: these should be included by developing a 
renormalization group treatment around mean-field the- 
ory, which is currently under construction [14 1 . 



To make our approach concrete, we study an assem- 
bly of N spherical particles of diameter a enclosed in a 
volume V in three spatial dimensions, interacting with a 
soft repulsion of finite range. To fix ideas we choose 



y(r <cr) = e(l-r/cr)", y(r > cr) = 0, 



(1) 



with r the interparticle distance, e the strength of the 
repulsion, and a = 2 (harmonic repulsion). Although 
several systems are described by a Hertzian repulsion 
(a = |), the harmonic model originally proposed to de- 



scribe wet foams jl5j has become a paradigm in numer- 
ical studies of the T = jamming transition 0, 0- It 
was also studied at finite temperatures [H, [3l i and finds 
experimental realizations in emulsions and soft colloids. 
The choice a = 2 is also technically more convenient, 
but we emphasize that our approach is easily generalized 
to any repulsive potential. The model has two control 
parameters: the temperature T, and the fraction of the 
volume occupied by the particles in the absence of over- 
lap: if — TiNa'^ /{dtV). We set tr and e to unity. 

Over the last decade, a number of numerical observa- 
tions were reported for this model A jamming tran- 
sition is observed at T = at some critical volume frac- 
tion, Lpj , the density above which packings carry a finite 
density of particle overlaps. Numerically, energy density, 
Cgs, and pressure, P, are found to increase continuously 
from zero above ipj as power laws The pair corre- 
lation function of density fluctuations gir), develops 
singularities near ipj [l^ , which are smoothed by thermal 
excitations [isj . In particular, g{l) = oo at ipj and T ~ Q, 
which implies that the density of contacts between parti- 
cles, z, jumps discontinuously from to a flnite value, Zr, 
at ipj . Above ipj , z increases algebraically with 0, Q ■ 
Thus, jamming appears as a phase transition taking place 
in the absence of thermal motion, with a peculiar critical 
behaviour and observable physical consequences 

The success of our approach relies on its ability to ac- 
curately describe dense systems of harmonic spheres at 
very low T, which is theoretically challenging [l7| ■ Simple 
liquid state theories, such as integral equations work 
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FIG. 1: Sketch of the derivation of the replicated free energy 
and effective potential in Eq. Q. Each particle in the original 
liquid is replicated m times (dashed spheres). Assuming that 
the replicated particles form a "molecule" of average cage size 
A, we trace out in the partition sum the degrees of freedom 
of (m — 1) copies of the liquid to obtain an effective one- 
component liquid (black spheres) interacting with an effective 
pair potential Kff(r) (green lines). 



well for dense systems only when T is not too low |18| . 
At lower T , numerical simulations [lij indicate the ap- 
pearance of a complex free energy landscape associated 
with slow dynamics, as found in glass-forming liquids. 
Theoretically, the mode-coupling theory of glasses can be 
applied [lij but only gives limited dynamical insights, in 
particular failing to identify the jamming transition |17| . 
For structure, a possible path is a recently developed an- 
alytical approach based on replica calculations [ll|, [2^ . 
While generic in principle, the method requires in prac- 
tice specific approximations. We have found that prior 
works on Lennard- Jones 11 1 and hard sphere Q models 
fail when directly applied to models of harmonic spheres 



near jamming 17|. The main technical achievement of 



the present paper is the derivation of a new analytical 
scheme to compute the structure and thermodynamics of 
the model ([T]) over a range of parameters broad enough 
to allow the study of the T = jamming transition. 

We introduce m copies of the system of harmonic 
spheres as a mathematical tool to probe its complex free 
energy landscape (20| . and develop approximations to 
study the statistical mechanics of the replicated liquid 
using an effective potential approach Q, as sketched in 
Fig. [TJ We make a Gaussian ansatz for the probability 
distribution of replicated particles around the center of 
mass of the "molecules" shown in Fig. [1] p(xi • • • Xm ) ^ 
/d3xn::Li(2^^)"^/'e-(""-^)'/(2A)^ which defines the 

cage size A. Our central approximation is now per- 
formed, in which only two-body interactions between par- 
ticles in copy 1 induced by the coupling to the other 
{m — 1) copies are retained, see Fig. [TJ Consider two 
"molecules" , each composed of m particles with positions 
(xi • • • Xm) and (yi • • • ym): the effective potential be- 



tween the particles of replica 1 is obtained by averaging 
the total interaction X^aLi ^(^a — Ya) over the positions 
of particles within the (m — 1) remaining replicas: 



-/3Ktf(xi-yi) 



X <^ p(xi • • •x,„)p(yi 



on 



-/3y(x„-y„) 



(2) 



Thanks to the Gaussian form of the integral the latter 
expression can be rewritten as follows: 



-/3Kff(r) _ 



du 



(3) 



where q{u) = / d^t e-/5^(«-*) e-*'/(4A)/(4^^)3/2 j^^g 
explicit expression in terms of error functions, and /3 = 
1 /T. Finally, the free energy F{m, A; ip, T) is obtained by 
considering K;ff (r) — mV{r) as small, which becomes ex- 
act when A — >• 0, see Eq. ([3]), and doing standard pertur- 
bation theory Q around the liquid with potential mV{r), 
which is equivalent to a liquid with potential V{r) at tem- 
perature T /m. We obtain an effective one-component 
system with a free energy parametrized by A and m: 



F{m,A;ip,T) 



T 

Fhm-m{m, A) + Fiiqif, — ) 
TO 



r (4) 

dr5iiq(r,^,^)[e-^[^-W-™^Ml-l], 

TO 



where ^iiq((p, T) and 5iiq(r, ip, T) are respectively the free 
energy and pair correlation function of the original (non- 
replicated) fluid, and i^harm = ~^[(to — l)ln(27ryl) + 
TO — 1 -|- In to] is the ideal gas contribution for the repli- 
cated system |ll| . Thus the core of the approximation 
is embodied by the effective potential Vcs (r) . Physically, 
the presence of (to — 1) replicas induces near jamming 
a strong short-range effective attraction, similar in spirit 
to depletion forces in colloid-polymer mixtures 

Our task becomes the study of a complicated effective 
fluid described by Eq. To simplify calculations we 
perform a standard approximation. 



5iiq(r, ^, T) EE e-''^Wy(r, ^, T) « e"^^Wy(l, ^, 0), (5) 

well-suited to study the T limit To compute 
-Fliq and j/(l, Lp, 0) analytically we choose the hypernetted 
chain approximation |3| , although more elaborate closure 
relations 0] could be used. This could change slightly 
the location of the transition, but not its nature or the 
scaling predictions we derive. Finally, to obtain concrete 
results for a given state point (ip,T), we minimize the 
free energy with respect to the cage size A and, in the 
glass phase, to the replica number to [iot . 

We first determine the location of the transition be- 
tween the fluid and glass phases, signaled by the appear- 
ance of a free energy minimum with to < 1 [ill . [20| . 
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FIG. 2: Theoretical phase diagram of soft repulsive spheres. 
The glass transition temperature Tk separates the liquid and 
glass phases with Tk ~ {(p—^pK)"^ near ^pK ~ 0.577. At T" = 0, 
the glass jams under compression across (pgcp ~ 0.633, above 
which no glass state with no particle overlap exists at T = 0. 
Thus, the ground state energy gga increases continuously from 



as 



if - fgcp) 



see Fig. [2l A finite temperature glass transition, Tk, 
emerges continuously from zero above fK ~ 0.577, as 
Tk ^ {(fi ~ fK)^- We obtain the full thermodynamic 
behaviour in the glass phase (notably energy, pressure, 
specific heat, glass fragility) which compares qualitatively 
well with numerical results [l6j . In particular, the ground 
state energy and pressure remain zero across ipK, show- 
ing that just above ipK, T = glasses are not jammed. 
In these glassy states, like in a hard sphere crystal, parti- 
cles can vibrate near well-defined (but random) positions, 
and the system is not jammed Q. 

We now concentrate on the T — >■ limit at large vol- 
ume fraction in the glass phase. We obtain the ground 
state energy shown in Fig. [21 As found in simulations , 
it grows continuously from zero above a critical pack- 
ing fraction, Cgs ~ (<^ — (^Sgcp)^, so that the pressure in- 
creases linearly, P ^ {(p — fgcp)- The 'glass close pack- 
ing' [1], tpgcp, represents in our calculations the largest 
density where T = glasses with no particle over- 
lap exist. Within the present approximation we obtain 
fgcp = 0.633 > (fK- Thus, from the sole knowledge of 
V{r) in Eq. ([1]), our theory predicts the existence and 
location of a jamming transition deep in the glass phase, 
and accounts for its critical nature. Since a large num- 
ber of metastable states exist in the glass phase, our 
approach also directly explains the strong protocol de- 
pendence of the critical jamming density fj observed in 
simulations 0, [2l| , which get arrested in nonequilibrium 
amorphous states, and thus jam at a critical packing frac- 
tion (fij < tpgcp. However, the results we obtain near <Pgcp 
below are found to hold for any metastable glass, and 
therefore also hold near any protocol-dependent (pj. 
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FIG. 3: Evolution of the maximum of the pair correlation 
function near contact with T and ip. While (?max diverges on 
both sides of the transition at T = as gmax ~ |(^ — (pgcp\~^ , 
this divergence becomes a smooth maximum at finite T near 
the transition whose pos ition shifts as Vt (dashed line), as 
observed numerically [l^ . and experimentally (Ts" 



We now turn to the calculation of the pair correlation 
function near fgcp- Within our approximation, g(r) is 
directly related to the effective potential: 



5(r) = e-^^--Wy(l,^,0), 



(6) 



and comes as a direct result of the free energy minimiza- 
tion. We concentrate on the physics of interparticle con- 
tacts and thus focus on distances r « 1, in the vicinity 
of the jamming transition (T l,(p ~ </5gcp). 

At T = 0, we find that g{r) develops a diverging peak 
near contact, which obeys the following scaling law: 



r-1 



(7) 



where J-"± (x) are asymmetric scaling functions which de- 
pend on the sign oi 6(p = ip — Vgcp a-nd can be com- 
puted analytically. In particular, \ogT+{x) cx —x^ and 
J^-{x) oc x~^ when x ^ 1. Note that T-{x) can be de- 
rived using the hard sphere potential Q . The scaling ([7]) 
means that the peak height, (?max, diverges as |(5</'l^^ 
both sides of the transition at T = 0, see Fig. [3j while 
its width vanishes as \5p\. This behaviour was found in 
simulations [l^ . In Fig. 2] we show that numerical results 
not only obey the scaling form in Eq. ([7]), but also that 
the asymmetric shape of the scaling functions compares 
extremely well with our theoretical predictions. 

At very low, but finite temperature, the peak diver- 
gence is smoothed by thermal fluctuations, which was 
the focus of a recent study [l3| . In Fig. [3] we show the 



predicted smooth evolution of ^max when Pgcp is crossed 
at finite T. A nonmonotonic evolution with density is 
obtained, as in experiments 13|, \2/^. The position of 
the maximum evolves with T with a scaling in perfect 
agreement with numerical work [l3l |. We go further and 
predict the thermal broadening of ^(r), sec Fig.S] We ob- 
tain nearly perfect agreement of theory with simulations 
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is general enough that it can be systenLaticahy improved 
and generahzed to various models, such that new, or 
more precise predictions could be made, hopefully foster- 
ing more numerical or experimental work. While Bcrnal 
saw packings as simplified models for atomic liquids, it 
is equally useful to consider packings as a special class of 
disordered ground states. 

We thank P. Chaudhuri for providing the configura- 
tions used in numerical verifications of theoretical pre- 
dictions in Fig. m H. Jacquin acknowledges financial 
support from Capital Fund Management (CFM) Foun- 
dation. L. Bcrthicr is partially funded by ANR Dynhet 
and Region Languedoc-Roussillon. 



FIG. 4: The pair correlation near jamming predicted by the- 
ory (full lines) and measured in numerical simulations (sym- 
bols). Left panels: scaling behaviour at T = above (top) 
and below (bottom) the jamming transition showing the con- 
vergence of the first peak near r = 1 to a delta function with 
asymmetric scaling functions on both sides of the transition. 
Right panels: the first peak of the pair correlation broadens 
when T increases at constant (p above (top: 5(p = 2.8 • 10~*) 
and below the transition (bottom 5(p = —3.5 • 10~*). To ease 
visualization, we show the evolution of g{r) / g^z.^, where gmax 
can be read from Fig. |3l 

over several decades of temperatures with no adjustable 
parameter. The full scaling of g{r) near contact, as a 
function of T and 5(p near the jamming transition, is the 
main new achievement of the present work. 

Finally we obtain the number of contacts per parti- 
cle by integration: z = 2Aip Jg drr^g{r). The diverging 
peak described by Eq. ^ gives a discontinuous jump of 
z from to = 6, the celebrated isostatic value, at (/Jgcp, 
as observed [HQ 111 and already derived in Q. Above the 
transition we find z—Zc oc i^p—ipgcp)'^ with 7 = 1. The ex- 
ponent is in quantitative disagreement with the observed 
7=^0- Indeed, this exponent has been related to the 
presence of fluctuations [5| that are presumably not well 
captured by our mean-field theory, indicating that more 
detailed calculations (possibly based on renormalization 
group [31) should be developed to predict g{r) over a 
broader range of intcrparticle distances. 

Our results show that the fully nonequilibrium problem 
of soft particle packings, relevant to understanding the 
mechanical properties of many soft materials, can be suc- 
cessfully address using equilibrium statistical mechanics 
tools. As an application we have developed a many-body 
theory of the jamming transition of soft repulsive spheres 
which satisfactorily derives, from first principles, the ex- 
istence and location of a jamming transition, and some of 
its peculiar critical behaviour, and makes new predictions 
for correlation functions of jammed states. Our approach 
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